#ifndef __CURVEFITTING__
#define __CURVEFITTING__

#include <iostream>
#include <vector>
#include "Spline.h"

using namespace std;

enum SplineType {
    ppform,
    bspline
};


template <int N, SplineType splinetype>
class CurveFitting {
private:
    int n;
    int dim;
    vector<double> length;
    vector<vector<double>> knots;// interpolation knots with dim dimensions
    vector<vector<double>> curve;// fitted points with dim dimensions
    
public:
    CurveFitting(vector<vector<double>> _knots) : knots(_knots) {
        n = knots.size();
        dim = knots[0].size();
        
        // fill length with the formula from the handout
        length.push_back(0);
        for (int i = 1; i < knots.size(); i++) {
            double res = 0;
            vector<double> &x1 = knots[i];
            vector<double> &x2 = knots[i-1];
            for (int j = 0; j < x1.size(); j++) {
                res = res + (x1[j] - x2[j]) * (x1[j] - x2[j]);
            }
            length.push_back(length[i - 1] + sqrt(res));
        }
        curve_fit();
    }

    void curve_fit() {
        double dl = 0.01;// step size when fitting the curve
        for (int i = 0; i < dim; i++) {
            // get the i-th dimension of the knots
            vector<double> x_dim;
            for (int j = 0; j < n; j++) {
                x_dim.push_back(knots[j][i]);
            }

            if (splinetype == ppform) {
                if (N == 1) {
                    ppformSpline<1> ppspline_1(length, x_dim);
                    vector<double> nodes;
                    for (double k = 0; k <= length[n - 1]; k += dl) {
                        nodes.push_back(ppspline_1(k));
                    }
                    curve.push_back(nodes);
                } else if (N == 3) {
                    ppformSpline<3, Natural> ppspline_3(length, x_dim);
                    vector<double> nodes;
                    for (double k = 0; k <= length[n - 1]; k += dl) {
                        nodes.push_back(ppspline_3(k));
                    }
                    curve.push_back(nodes);
                } else {
                    cerr << "N should be 1 or 3" << endl;
                    exit(-1);
                }
            } else if (splinetype == bspline) {
                if (N == 1) {
                    Bspline<1> bspline_1(length, x_dim);
                    vector<double> nodes;
                    for (double k = 0; k <= length[n - 1]; k += dl) {
                        nodes.push_back(bspline_1(k));
                    }
                    curve.push_back(nodes);
                } else if (N == 3) {
                    Bspline<3, Natural> bspline_3(length, x_dim);
                    vector<double> nodes;
                    for (double k = 0; k <= length[n - 1]; k += dl) {
                        nodes.push_back(bspline_3(k));
                    }
                    curve.push_back(nodes);
                } else {
                    cerr << "N should be 1 or 3" << endl;
                    exit(-1);
                }
            }
        }
    }

    vector<vector<double>> get_curve() { 
        return curve; 
    }
};

#endif
